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Abstract 

A  time-reversal  mirror  is,  roughly  speaking,  a  device  which  is  capable  of 
receiving  an  acoustic  signal  in  time,  keeping  it  in  memory  and  sending  it  back 
into  the  medium  in  the  reversed  direction  of  time.  In  this  paper,  we  employ 
an  accurate  numerical  method  for  siinulating  waves  propagating  in  complex 
one-dimensional  media.  We  use  numerical  simulations  to  reproduce  the  time- 
reversal  self-averaging  effect  which  takes  place  in  randomly  layered  media. 
This  is  done  in  the  regime  where  the  inhomogeneities  are  smaller  than  the 
pulse,  which  propagates  over  long  distances  compared  to  its  width.  We  show 
numerical  evidence  for  possible  use  of  an  expanding  window  time-reversal 
technique  for  detecting  anomalies  buried  in  the  medium. 


1.  Introduction 

A  time-reversal  mirror  is,  roughly  speaking,  a  device  which  is  capable  of  receiving  an  acoustic 
signal  in  time,  keeping  it  in  memory  and  sending  it  back  into  the  medium  in  the  reversed 
direction  of  time.  In  the  context  of  ultrasounds,  time-reversal  mirrors  have  been  developed 
and  their  effect  studied  experimentally  by  Mathias  Fink  and  his  collaborators  at  the  Laboratoire 
Ondes  et  Acoustique  (ESPCI-Paris)  [8,  9].  The  main  effect  is  the  refocusing  of  the  scattered 
signal  after  time  reversal  in  a  random  medium:  an  acoustic  pulse  is  sent  in  a  disordered 
medium  generating  a  ‘noisy’  reflected  signal  which  is  time  reversed  and  sent  back  into  the 
medium.  The  new  reflected  signal  is  a  pure  pulse  with  a  shape  similar  to  the  initial  pulse. 
Amazingly,  its  ‘refocusing’  takes  place  in  time  and  space  and  seems  to  be  independent  of  the 
realization  of  the  medium,  in  the  regime  where  its  correlation  length  is  smaller  than  the  typical 
wavelength  of  the  pulse. 

From  the  theoretical  point  of  view,  a  first  proof  of  this  refocusing  effect  has  been  obtained 
in  [7]  in  the  context  of  a  one-dimensional  random  medium  for  which  only  the  time  refocusing 
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is  relevant.  The  refocusing  is  obtained  by  using  asymptotics  in  the  regime  where  there  are 
three  well-separated  scales: 

correlation  length  of  the  medium  ^  wavelength  of  the  pulse  ^  ^ 
wavelength  of  the  pulse  distance  of  propagation 

The  fluctuations  of  the  medium  are  not  assumed  to  be  small  but  rather  of  the  order  of  several 
tens  of  per  cent.  In  the  case  of  a  stationary  random  medium,  a  theoretical  formula  for  the 
limiting  refocused  pulse  is  given  in  [7].  It  shows  that  the  refocused  pulse  is  a  convolution  of 
the  initial  pulse  with  a  kernel  that  depends  on  the  cut-off  of  the  reflected  signal  and  on  the 
autocorrelation  of  the  medium.  In  the  non-stationary  case,  which  is  the  case  for  the  detection 
problem  considered  in  this  paper,  there  is  no  explicit  formula  available. 

In  this  study,  we  write  a  numerical  method  in  the  frequency  domain  for  the  one¬ 
dimensional  (ID)  wave  equation  and  use  it  to  study  the  time-reversal  refocusing  effect  for 
waves  propagating  through  a  scatterer  comprising  many  layers  with  randomly  varying  wave 
speed.  On  each  layer  of  the  scatterer,  our  method  yields  two  relations  in  terms  of  the  unknown 
field  and  its  normal  derivative  on  the  interfaces  of  that  layer.  When  matched  to  an  incident 
field  on  the  first  interface,  our  approach  is  analogous  to  a  transfer  matrix  method  [1]  that 
is  formulated  in  the  frequency  domain.  We  derive  our  ID  method  via  a  boundary  integral 
equation  approach  in  the  frequency  domain  to  facilitate  future  extension  to  higher  dimensions. 
For  example,  in  2D,  a  boundary  integral  method  has  been  developed  and  used  to  study  resonant 
modes  for  electromagnetic  scattering  in  photonic  crystal  slabs  [11].  A  frequency  domain 
boundary  integral  formulation  is  also  the  appropriate  context  for  implementation  of  a  fast 
multipole  method  [13].  Using  our  ID  method,  we  simulate  the  time-reversal  refocusing  effect 
in  reflection,  and  analyse  the  potential  for  using  time  reversal  in  detecting  the  presence  of  an 
extended  strong  scatterer  buried  in  the  random  medium.  The  inverse  problem,  which  consists 
of  reconstructing  the  large  scale  variations  of  the  medium  (its  slowly  varying  component),  is 
an  extremely  challenging  problem.  It  has  been  addressed  in  [1,  2]  and  [12]  by  using  reflected 
signals.  Theoretically,  the  slowly  varying  components  of  the  medium  parameters  appear  in 
the  coefficients  of  an  infinite  system  of  transport  equations  [1]  and  this  makes  the  inverse 
problem  very  delicate.  In  this  paper,  we  propose  to  tackle  numerically  a  simpler  problem 
(i.e.  detection  of  an  anomaly  of  size  larger  than  the  typical  inhomogeneities),  by  using  an 
expanding  window  time-reversal  technique  (section  4.2)  and,  in  particular,  by  exploiting  the 
self-averaging  property  of  the  refocused  pulses. 

Time  reversal  in  several  dimensions  has  been  studied  in  various  regimes:  for  instance,  in 
[5]  in  the  parabolic  approximation  regime,  in  [10]  in  the  case  of  a  point  source  in  a  randomly 
layered  medium  and  in  [3]  for  general  waves  propagating  in  weakly  fluctuating  media.  Recent 
applications  to  imaging  can  be  found  in  [6].  The  numerical  method  presented  here  in  the 
context  of  strongly  contrasted  media  is  well  adapted  to  multi-dimensional  generalizations 
which  are  part  of  our  future  work. 


2.  Formulation  of  interface  equations  in  the  frequency  domain 

Our  numerical  formulation  is  based  on  a  frequency  domain  boundary  integral  method  which, 
in  the  ID  context,  yields  a  transfer  matrix  for  the  unknown  components  of  the  field  and 
its  normal  derivative  on  each  layer  interface.  While  in  ID,  our  method  is  not  particularly 
novel,  the  frequency  domain  boundary  integral  formulation  can  be  extended  to  simulate  time 
reversal  in  higher  dimensions,  where  theory  is  available  only  in  special  cases.  Extension  of 
the  numerical  method  is  briefly  discussed  in  section  2.5. 
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21.  Green’s  functions  in  the  fi‘equency  domain 


We  consider  the  one-dimensional  wave  equation  for  a  field  u{x,t)\ 


d^u 


X  e  Q,  t  >  0 


(1) 


where  Q  =  (— oo,  oo)  and  c{x)  is  the  wave  speed.  Let  Q  —  (— oo,  cq)  L  Qq  LI  •  ■  •  U  LJ 
{On,  oo)  where  Qj  ==  (aj,  aj+i),  j  =  0, . . . ,  1  (ao  <  ai  <  •  •  •  <  on).  A  layered  medium 

is  considered  by  taking  c(x)  to  be  piecewise  constant  on  Q  as  follows: 


fcj  forjc  e  Q;,  j  =0,  1  _ 

~  \c'  foTx  e  (-00,  ao)  U  (^a,,  oo).  ^  ^ 

We  consider  the  case  in  which  the  continuous  quantities  on  the  interfaces  x  ~  uq,  a 
are  u  and  Denote  by  P(co)  tlie  Fourier  transform  in  time  of  a  function  /(/)  (Fico)  = 

/_“  /(Oe^'^dr). 

Taking  the  Fourier  transform  of  (1)  using  (2),  we  obtain  the  following  Helmholtz  equations 
for  each  frequency  a>  (—oo  <  co  <  oo): 


X  e  (kj^^co/cj) 

(3) 

x  €  (-OC,  ao)  U  (a^v,  oo)(k  =  w/c'). 

(4) 

Introduce  a  Green’s  function  H’'  (x,  y)  =  (2ikj)~‘^  for  each  layer  ilj  of  the  medium, 

and  a  free  space  Green’s  function  H(x,  y)  =  (2ik)~^  Qik\x-y\^  These  Green’s  functions  give 
rise  to  outward  radiating  waves  at  jc  =  -oo  and  satisfy  equations  (3)  and  (4),  respectively, 
with  a  non-homogeneous  singular  source  term  ^  (x  —  y),  where  5  (jc)  is  the  Dirac  delta  function. 


2.2.  Interior  representations 


To  represent  the  field  u  in  the  interior  of  layer  Qy,  we  apply  reciprocity  to  (3),  integrating  on 
the  interval  Q,j  —  {aj,aj+\)  to  obtain,  for  y  e  (aj,  ay+i), 

Using  the  filtering  property  of  5(jc)  and  integration  by  parts,  the  domain  integrals  cancel  out 
and  (5)  reduces  to 


u{y)  = - - H(ay+i) - - u(aj)  -  y) — — — 


+  Nj(aj,y) 


du(aj) 

3^’ 


y  €  (aj,aj+i). 


(6) 


Equation  (6)  determines  the  field  u  in  the  interior  of  layer  Q j  given  the  (as  yet  unknown) 
values  of  the  field  w  and  its  normal  derivative  |”  at  the  boundaries  of  the  layer. 

To  obtain  an  interior  representation  to  the  left  of  the  scatterer,  consider  (6)  on  the  interval 
(Oj,  aj+i)  =  (— L,  ao)  using  the  Green’s  function  H(x,y),  As  our  incident  wave,  we  choose 
a  plane  wave  propagating  from  left  to  right.  Hence,  the  field  to  the  left  of  the  scatterer  can  be 
written  as  w  =  e'^^'  +  Wsc,  where  Msc  is  the  scattered  wave  on  (~L,  ao).  Substituting  u  into  (6), 
and  using  the  relations  H(-L,y)  =  (2i^)"*  and  the  second 

and  fourth  terms  of  (6)  reduce  to  e**'' .  Taking  the  limit  as  L  oo,  an  interior  representation  for 
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the  reflected  field  to  the  left  of  the  scatterer  is  obtained: 

uiy)  =  -  H(ao,  y  6  (-L,  ao).  (7) 

OX  oy 

Similarly,  by  considering  (6)  on  the  interval  L),  where  there  is  no  incident  wave,  the  limit 
as  L  00  yields  the  following  interior  representation  for  the  transmitted  field  to  the  right  of 
the  scatterer: 


^(y)  = - - u(a!^)  +  H{aN,  y)- 


Sx 


By 


y  e  (aN^oo). 


(8) 


23.  Formulation  of  a  system  of  equations  in  the  interface  unknowns 

It  remains  to  determine  the  unknown  values  of  the  field  it  and  its  noimal  derivative  ~  on  the 
interfaces  y  =  ao^ai, . . .  Using  a  boundary  integral  method,  we  obtain  a  closed  system 
of  linear  algebraic  equations  for  these  unknowns  by  taking  the  limit  as  the  interior  point  y 
approaches  the  two  boundary  points  of  each  layer. 

Consider  the  interior  representation  (6)  in  layer  of  the  medium.  Taking  the  limit  as 

y  -V  at,  and  using  the  fact  that  =  -1/2  and  we  obtain  a 

relation  in  terms  of  the  interface  unknowns  for  layer  (j  =  0, 1 , . . . ,  ^  1) : 


1  du(aj) 
2ikj  By 


Bx 


uioj+i)  -  H^{aj+ua'^) 


Bu{ttj+\) 

By 


(9) 


Similarly,  taking  the  limit  of  (6)  as  y  and  using  the  fact  that  — ^  zz:  1/2  and 

Hl(aj+\ ,  we  obtain  a  second  relation  involving  the  interface  unknowns  for 

layer  (j  =  0, 1, . . . ,  N  —  1): 


2ik 


Bu(aj+^)  _ 

ay 


BHl{aj,aj^^)  _  ^Bu{aj) 


Bx 


(10) 


Equations  (9)  and  (10)  comprise  2N  equations  in  the  2Ar  +  2  interface  unknowns.  The 
r^naining  two  equations  are  obtained  by  perfonning  a  similar  limiting  process  on  intervals  to 
the  left  and  right  of  the  layered  medium  to  obtain,  respectively: 


-Hoo)- 


1  aM(flo) 

m  By 


and 


-M(ajv) 


1  Bu(aj^) 
m  By 


=  0. 


(11) 


The  closed  system  of  equations  for  the  2iV  +  2  unknown  values  of  the  u  and  on  the 
interfaces  located  at  y  =  ao,  ...,x  =  consists  of  (9),  (10)  (for  j  =  0,  1, . . . ,  AT  —  1) 
and  (11).  This  system  of  equations  is  assembled  into  a  banded  matrix  (bandwidth  5)  that  is 
effectively  a  transfer  matrix  for  the  unknown  components  of  the  field  and  its  normal  derivative 
on  the  layer  interfaces.  Solutions  of  the  linear  system  are  obtained  at  each  frequency  o)  using 
a  banded  solver  which,  for  N  layers,  requires  O  (N)  operations. 


2.4.  Example:  calculating  reflection  and  transmission  amplitudes 


Once  the  interface  unknowns  are  determined  via  solution  of  (9)— (11)  at  a  particular  reduced 
frequency  k  =  &)/c',  the  transmitted  wave  on  the  interval  y  €  oo)  is  determined  from  (8) 
as 


u^iy)  =  f{k)e^y. 


where 


t{k)  = 


uia^)  + 


1  ‘du(aN)\ 
ik  dy  )' 


(12) 
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Similarly,  the  reflected  wave  on  the  interval  y  €  (-00,  ^o)  is  determined  from  (7)  as 
Uzcfiy)  =  Mk)  +  e“•^  where  Rik)  =  —  (^m(£Jo)  -  —  j  .  (13) 


2.5.  Extension  of  the  boundary  integral  method  to  higher  dimensions 

While  in  ID,  equations  (9)-(ll)  yield  a  banded  transfer  matrix,  their  formulation  using  a 
boundary  integral  approach  extends  to  the  development  of  numerical  methods  for  multiple 
scattering  in  higher  dimensions.  For  example,  in  higher  dimensions  the  layers  Q.j  can  be 
closed  regions  that  define  scatterers,  with  prescribed  geometric  and  material  properties,  that  are 
embedded  in  a  homogeneous  background  medium.  In  tlie  interior  representation  (6),  the  right- 
hand  side  becomes  a  line  integral  (2D)  or  surface  integral  (3D)  written  in  terms  of  the  field  and 
its  normal  derivative  on  the  scatterer  interfaces.  The  ID  Green’s  functions  for  the  Helmholtz 
equation  are  replaced  by  higher  dimensional  counterparts  interior  and  exterior  to  the  scatterers. 
The  limiting  process  of  section  2.3  yields  a  closed  system  of  boundary  integral  equations  on 
either  a  collection  of  closed  contours  (2D)  or  closed  surfaces  (3D).  These  integral  equations 
are  solved  numerically  using  a  collocation  method  to  obtain  approximations  to  the  field  and 
its  normal  derivative  on  the  scatterer  interfaces.  Equation  (12)  for  the  transmitted  field  and 
equation  (13)  for  the  reflected  field  then  involve  integrals  written  in  terms  of  the  computed 
field  approximations  over  the  union  of  scatterer  interfaces.  In  higher  dimensions,  boundary 
integral  methods  reduce  the  computational  mesh  by  one  dimension  at  the  expense  of  a  fully 
populated  linear  system  in  the  collocation.  However,  analytical  manipulation  of  the  Green’s 
functions  can  be  used  to  uncouple  near  and  far-field  effects,  leading  to  sparse  linear  systems 
and  significant  acceleration  of  the  numerical  solution  of  the  governing  boundary  integral 
equations. 

Using  the  approach  described  above,  2D  boundary  integral  numerical  methods  were 
developed  for  analysis  of  bandgaps  and  resonances  in  2D  photonic  crystal  slabs  (finite  in 
X,  periodic  in  F,  homogeneous  in  Z)  [4,  11,  14].  These  methods  were  used  to  conduct  a 
parametric  analysis  of  resonance  quality  for  Fabry-Perot  laser  cavities  comprising  circular 
digjieciric  rods  in  both  the  lossless  [14]  and  lossy  [4]  cases.  Via  analytical  manipulation  of 
the  2D  Green’s  function,  the  formulation  was  accelerated  in  [11]  to  give  a  method  that  was 
0{N)  in  the  direction  of  finite  extent  (X),  where  N  is  the  number  of  interface  collocation 
points.  This  fast  method  was  used  to  compute  complete  photonic  bandgaps  (i.e.  for  ail  incident 
angles)  with  random  perturbations  in  crystal  geometry,  and  to  analyse  resonance  quality  in  the 
presence  of  channel  defects.  Formulations  for  the  development  of  efficient  boundary  integral 
numerical  methods  in  the  general  2D  and  3D  cases  are  presented  by  Rokhlin  in  [13].  These 
studies  can  serve  as  the  basis  for  numerical  simulation  of  time  reversal  in  higher  dimensions, 
where  theoretical  results  are  available  only  in  special  cases.  However,  efficient  and  accurate 
analysis  in  higher  dimensions  poses  a  challenging  computational  problem  due  to  the  inherent 
random  and  multi-scale  nature  of  the  time-reversal  effect. 

3.  Time  reversal  in  the  frequency  domain 

We  now  formulate  time  reversal  of  a  reflected  pulse  in  terms  of  the  reflection  scattering 
operator  R{k)  given  in  (13).  Consider  a  unit-amplitude  Gaussian  derivative  pulse  WincCy,  0  = 
y-(>’-yo-cf)  centred  at  y  =  yo  when  t  =  0,  where  f(x)  =  ~^x€  e*^^  and  e  is  a 

characteristic  wavelength  of  the  pulse.  With  this  normalization  the  maximum  value  of  |/|  is 
exactly  1 .  The  Fourier  transform  of  the  pulse  is 


190 


M  A  Haider  et  al 


Uinciy,  CO)  =  -  (14) 

c' 

=  A(a)/c')e“^/‘^',  where  i(fe)  =  v^e^ifee-**^<>e-^'*'/^  (15) 

By  (13),  the  reflected  signal  to  the  left  of  the  scatterer  is  given  by 

UrAy.  oi)  =  A{cold){R{w/c')  e-“>’/^'  +  e*“^^"').  (16) 

We  now  sample  a  portion  of  the  reflected  signal  using  a  window  function  Gi^{t)  — 
H{t)  -  H{t  -  to),  where  H{t)  is  the  Heaviside  step  function.  The  window  function  is 
used  to  record  a  sample  of  the  reflected  wave  at  y  =  y*  starting  at  time  t  =  t*.  Hence,  the 
reflected  wave  is  sampled  for  a  duration  to  in  the  time  interval  t  €  (r*,  r*  +  ?o)-  The  sample  of 
the  reflected  signal  in  the  time-domain  Wref  (0  is  given  by 

MrefCO  =  Mref(y,  f {y  -  y*)/cO  (17) 

where  the  Fourier  transform  of  Mref(y,  0  is  given  by  (16).  The  sampled  signal  (17)  is  then 
reversed  in  time  to  obtain  a  new  incident  wave  Mtr(0  =  “ref  (^*  +  Using  (17),  the  new 

incident  wave  Mo-evCy,  0  can  be  written  in  the  frequency  domain  as  the  convolution: 

fiuev(y,  CO)  = 

=K[e“('‘-"'“'">«/^)(-iw)N/2^e^  +  eT^y/'^)]  (18) 

where  R{a)/c')  is  the  complex  conjugate  of  R(o)/d).  By  (13),  the  second  reflection  of  the  new 
(time-reversed)  incident  wave  is  then  obtained  as  (o)  =  «trev(>S  co){R{o)/c')  + 

which  can  be  written  explicitly  as 

()-,  CO)  =  —^[^{R(co/c')  +  e*"5'A')  f”  cb 

c'  y  27V  J-oo 

p}(o>—Q})to  _  1 

i(a)  —  co) 

Fourier  inversion  of  (19)  gives  the  time-reversed  pulse  t)  in  the  time  domain. 

4>NiimericaI  simulations 

We  employed  the  method  described  in  section  2.3  in  the  numerical  simulation  of  time  reversal 
for  the  wave  equation  (1)  and  (2)  with  wave  speeds  cj  in  the  layers  varying  randomly 
about  a  background  velocity  c'  =  1.  The  layered  random  medium  was  located  on  the 
interval  [0, 1]  and  layers  of  a  fixed  width  l/N  were  considered  by  placing  the  interfaces  at 
aj  =z  j/N  ij  =  0, 1, . . . ,  N).  The  initial  pulse  location  yo,  signal  recording  location  y  and 
time-reversal  recording  location  y*  were  all  taken  to  be  one  unit  to  the  left  of  the  layered 
medium  (i.e.  yg  =  y  =  y*  =  *-1).  Since  c  =  c'  =  1  on  the  left  of  the  layered  medium, 
the  time  at  which  sampling  of  the  reflected  signal  commences  was  taken  as  t*  =  2.  The 
remaining  parameters,  which  control  the  multiple  scales  in  the  problem,  are  the  characteristic 
pulse  wavelength  e,  the  number  of  layers  N  and  the  length  of  the  sampling  interval  /q-  This 
set-up,  used  in  the  detection  problem  in  section  4.2,  corresponds  to  an  expanding  time  window 
with  t*  fixed  and  to  increasing.  Alternatively,  one  can  consider  a  sliding  time  window  which 
would  correspond  to  t*  increasing  and  to  fixed. 

For  a  particular  random  realization  of  the  layered  medium  we  first  computed,  at  each 
frequency,  the  reflection  amplitude  by  solving  the  linear  system  (9)-(l  1)  and  using  the  result 
in  (13),  We  considered  a  finite  number  of  equally  spaced  frequencies  and,  by  symmetry,  it  was 
sufficient  to  evaluate  all  functions  in  the  frequency  domain  at  only  non-negative  frequencies. 
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Figure  1.  A  pulse  scattered  by  a  randomly  layered  medium  generates  a  noisy  reflected  signal 
shown  in  (a),  and  magnified  in  (b).  A  simulation  with  five  random  realizations  (c)  shows  that  the 
reflected  signal  depends  on  the  realization  of  the  medium.  The  reflected  signal  is  sampled  using  a 
window  of  size  /q,  time  reversed  and  sent  back  into  the  same  medium.  Time-reversal  refocusing 
is  demonstrated  in  (d)  for  five  realizations,  as  well  as  the  self-averaging  property  of  the  refocused 
pulse. 


All  computations  were  performed  using  a  1  GHz  dual-processor  Apple  G4  Power  Macintosh 
in  MATLAB  6.5  which  employs  16-digit  precision.  All  simulations  were  conducted  for 
N  =  1000  layers  and  a  banded  solver  was  used  for  rapid  and  accurate  inversion  of  the 
coefficient  matrix  in  the  assembled  linear  system.  The  time-reversed  pulse  was  calculated 
in  the  time  domain  via  numerical  integration  of  the  convolution  integral  in  (19)  followed 
by  a  second  numerical  integration  for  Fourier  inversion.  In  both  cases,  Simpson’s  rule 
was  employed.  The  same  set  of  sampling  frequencies  were  used  in  computation  of  the 
scattering  operator,  and  in  the  numerical  integration  of  (19)  and  the  Fourier  inversion  formula. 
Once  the  scales  for  the  layers  and  pulse  were  set  via  N  and  6,  respectively,  the  frequency 
resolution  required  for  numerical  convergence  could  be  determined.  In  the  finest  scale  regime 
considered  in  this  study  (i.e.  N  =  1000,  6  =  O.l/ViV),  we  found  that  the  frequency  sampling 
CO  e  (-1000, 1000),  Aco  =  0.25  was  sufficient  for  visual  overlapping  of  successively  refined 
computations  of  the  time-reversed  pulse  in  the  time  domain. 


4,L  Pulse  refocusing  and  self-averaging  via  time  reversal 

Numerical  simulation  of  time  reversal  in  reflection  is  illustrated  in  figure  1  for  the  case 
N  =  1000  with  €  —  r  =  10  and  Cj  ==  1/^1  +  where  rjj  is  o.  sequence  of  i.i.d. 


Figure  2.  Refocusing  and  self-averaging  of  the  time-reversed  pulse  is  demonstrated  as 
the  characteristic  pulse  wavelength  e  =  is  varied  via  the  parameter  r  in  the  case  N  =  1000. 

The  factor  r  compensates  for  a  small  autocorrelation  produced  by  our  i.i.d.  medium.  It  is  fixed  in 
the  regime  where  the  multiple  scales  become  well  separated,  and  from  our  experiments  (a)-{d)  the 
value  r  =  10  shows  a  good  tuning  of  the  pulse  to  the  medium  in  terms  of  quality  of  self-averaging 
of  the  refocused  pulse.  For  each  value  of  r,  the  refocused  pulse  is  shown  for  ten  random  realizations 
of  the  layered  medium. 

random  variables  that  are  uniformly  distributed  over  the  interval  (—0.8,  0.8).  Observe  that 
in  the  homogenization  regime  the  quantity  to  be  homogenized  is  l/c^{x),  so  that  the  way 
we  modelled  the  velocity  random  fluctuations  corresponds  to  centred  fluctuations  around  the 
homogenized  quantity  1  /c^  =  1.  The  time  required  to  compute  the  refocused  pulse  in  the 
time  domain  for  a  single  realization  of  the  layered  medium  was  roughly  14  min. 

At  the  recording  location  y  =  —  1,  a  typical  time  signal  demonstrating  reflection  of  the 
incident  pulse  is  shown  in  figure  1(a).  Using  the  window  function  G/^ ,  a  sample  of  the  reflected 
signal  of  duration  to  is  recorded  (figure  1(b)),  reversed  in  time  and  then  re-introduced  into  the 
layered  medium  as  a  new  incident  wave.  Multiple  realizations  demonstrate  that  the  reflected 
wave  does  not  exhibit  a  self-averaging  property  that  can  capture  characteristics  of  the  layered 
medium  (figure  1(c)).  In  contrast,  the  reflection  of  the  new  (time-reversed)  incident  wave 
refocuses  the  pulse  in  a  manner  that  is  self-averaging  witli  respect  to  random  realizations  of 
wave  speed  in  the  component  layers  of  the  scattering  medium  (figure  1(d)). 

Tn  figure  2,  we  illustrate  self-averaging  and  refocusing  of  the  time-reversed  pulse  as  its 
characteristic  wavelength  e  =  is  tuned  to  the  fine  scale  of  the  random  medium  via  the 
parameter  r.  For  each  r,  we  computed  the  time-reversed  pulse  for  ten  random  realizations  of 
the  layered  medium.  As  the  pulse  wavelength  is  decreased,  separation  between  the  fine 
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scale  of  random  variation,  the  intermediate  scale  of  pulse  wavelength  and  the  distance 
of  propagation  leads  to  exhibition  of  the  self-averaging  property  of  the  refocused  pulse 
when  r  =  10  (figure  2(d)).  In  this  regime,  with  to  relatively  large  so  that  the  recording 
window  contains  most  of  the  reflected  signal,  the  physical  wavelength  of  the  refocused  pulse 
matches  that  of  the  original  incident  pulse  with  roughly  a  12%  loss  in  pulse  amplitude  due 
to  transmission  of  waves  through  the  layered  medium  not  yet  reflected.  Our  finding  is  in 
agreement  with  the  theoretical  formula  obtained  in  [7]  which  shows  that,  asymptotically  in 
the  regime  of  separation  of  scales  (e  0),  the  refocused  pulse  ftr  is  obtained  as  a  convolution 
f^(t)  =  (/(■)  ★  K{‘))  (t)  of  the  time-reversed  initial  pulse  /  with  a  kernel  K  given  by  its 
Fourier  transform  K(a))  =  A(o),  r)  dr.  The  density  A  is  known  explicitly  in  die  case 

of  homogenous  media.  We  hypotliesize  that  this  self-averaging  property  of  the  time-reversed 
refocused  pulse  can  be  employed  in  detecting  the  presence  of  a  strong  homogeneous  scatterer 
that  is  buried  in  the  layered  random  medium.  The  idea  is  that  to/2  corresponds  to  a  depth 
of  penetration  and  therefore  a  comparison  of  the  kernels  K  with  and  without  a  buried  object 
produces  a  detection  tool.  We  choose  here  not  to  rely  on  formulae  or  computation  of  the 
densities  A,  since  these  will  not  be  available  in  the  multi-dimensional  context.  Rather,  we 
use  a  comparison  between  maxima  (amplitudes)  of  the  refocused  pulses  with  and  without  the 
buried  object. 

For  practical  applications,  we  suggest  that  the  experiment  be  repeated  at  various  locations 
in  order  to  acquire  maxima  of  the  refocused  pulse  without  the  buried  object  as  a  function  of  to, 
and  dicn  test  other  locations  by  comparison  as  described  above.  Strictly  speaking  this  is  not 
compatible  with  the  assumed  layered  structure  but  one  can  consider  that  for  remote  locations 
the  random  layering  is  independent. 


4.2.  Detection  of  a  buried  extended  strong  scatterer  via  time  reversal 

Numerical  simulations  were  conducted  in  the  self-averaging  regime  described  in  the  previous 
section.  We  hypothesized  that  time-reversal  detection  can  be  evaluated  by  comparing  the 
distribution  of  refocused  pulse  amplitudes  over  multiple  random  realizations  of  the  layered 
medium  in  the  presence  and  absence  of  a  buried  strong  scatterer.  The  strong  scatterer  is 
introduced  by  setting  the  wave  speed  for  L  layers  at  a  location  y  =  to  50c\  where  d  is  the 
wave  speed  in  the  homogeneous  background  medium.  The  case  L  —  100,  yc  =  0.67  is 
illustrated  in  figure  3(a).  A  sample  realization  of  the  reflected  signal  for  the  time-reversal 
parameters  t*  =  2.0,  to  =  3.0  is  shown  in  figure  3(b).  We  observe  that  the  reflected  signal 
contains  no  clear  signature  indicating  the  presence  of  a  strong  scatterer  buried  in  the  layered 
random  medium. 


4 .2 .7 .  Detecting  the  presence  of  the  scatterer.  We  analysed  the  effect  of  the  strong  scatterer  on 
the  amplitude  of  the  time-reversed  refocused  pulse  for  multiple  time-reversal  recording  times 
to  (figure  3(c)).  An  expanding  window  time-reversal  analysis  was  conducted  by  increasing 
the  parameter  to  increments  of  0.5  s  up  to  a  value  of  3.0  s.  For  each  recording  time, 
we  considered  ten  random  realizations  with  and  without  a  buried  strong  scatterer.  All  100 
random  realizations  in  figure  3(c)  were  independent.  As  the  recording  time  to  is  increased  to 
incorporate  portions  of  the  reflected  signal  that  correspond  to  the  depth  of  the  buried  object 
(i.e.  to  =  1.5  s),  we  observe  a  statistically  significant  difference  in  mean  pulse  amplitude 
between  the  cases  with  and  without  a  buried  strong  scatterer.  The  direct  coherent  reflection 
from  the  strong  scatterer  interface  (where  c  increases  abruptly  from  1  to  50)  is  responsible 
for  this  enhancement  of  the  amplitude  of  the  refocused  signal.  However,  due  to  the  strong 
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Random  Layered  Medium 


N=1000  layers 


Random  Layered  Medium  with  Strong  Scatterer 


(a) 


c=c'=1  c=50c’ 


c=:c'=1  N=1000  layers 


Reflected  signal 
(no  strong  scatterer) 


Figure  3.  Numerical  simulations  for  detection  using  an  expanding  window  time-reversal  analysis 
{N  —  1000  layers,  €  =  (a)  Independent  realizations  of  the  medium  with  and  without  a 

buried  strong  scatterer.  A  pulse  is  incoming  from  the  right  and  the  scatterer  consists  of  L  =  100 
layers  at  a  depth  yc  =  0.67  with  wave  speed  50c'.  (b)  Reflected  signals:  the  coherent  reflection 
due  to  the  strong  scatterer  is  hidden  in  the  noise  due  to  the  inhomogeneities,  (c)  Amplitudes  of 
refocused  pulses  as  the  time-reversal  window  expands  (5  windows,  20  independent  realizations 
per  window:  10  witliout  buried  object  and  10  with  buried  object).  Horizontal  lines  indicate  mean 
amplitudes  for  each  subset  of  ten  realizations. 

fluctuations  in  the  mediuni,  this  coherent  reflection  was  hidden  in  the  reflected  signal,  and,  in 
that  sense,  time  reversal  reveals  it. 

Self-averaging  of  the  pulse  amplitude  improves  as  the  signal  recording  time  is  increased 
to  fo  —  3.0  s.  In  this  self-averaging  regime,  these  results  indicate  the  potential  use  of  the 
time-reversed  refocused  pulse  as  a  measure  for  detection  of  a  buned  strong  scatterer. 

To  assess  the  accuracy  of  the  proposed  detection  measure,  we  computed  the  distribution 
(i.e.  mean,  standard  deviation)  of  the  refocused  pulse  amplitude  in  the  case  to  —  3.0  s  for  50 
random  realizations  with  and  without  the  presence  of  the  strong  scatterer  described  above. 
The  effect  of  scatterer  size  was  quantified  by  varying  the  parameter  L  and  the  random  number 
generator  was  assigned  a  new  seed  for  each  value  of  L.  The  computed  probability  distribution 
functions  are  shown  in  figure  4. 

For  the  case  of  no  scatterer  (L  =  0) ,  the  mean  pulse  amplitude  was  0.8786  with  a  standard 
deviation  of  0.0132  while  for  the  largest  buried  object  considered  (L  =  100)  the  mean  pulse 
amplitude  was  0.9396  with  a  standard  deviation  of  0.0151. 

For  instance,  by  choosing  a  detection  amplitude  threshold  at  the  level  0.91,  we  find  that: 

•  the  probability  of  error  hy  false  alarm,  that  is  the  area  to  the  right  of  0.91  imder  the  L  =  0 

probability  density,  is  less  than  1%; 
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Figure  4.  Probability  distributions  of  refocused  pulse  amplitude  with  varying  object  size  L. 
For  each  buried  strong  scatterer  of  size  L,  the  refocused  pulse  amplitude  was  computed  for  50 
independent  random  realizations  of  the  layered  medium.  The  mean  and  standard  deviation  of 
the  amplitudes  were  fit  to  a  nomial  distribution,  (a)  Wlien  the  buried  strong  scatterer  is  as  small 
as  5%  of  the  size  of  the  layered  medium  {L  —  50),  the  probability  of  error  in  the  time-reversal 
detection  measure  is  small,  (b)  When  the  size  of  the  buried  strong  scatterer  is  below  5%  of  the  size 
of  the  layered  medium,  the  probability  of  error  in  the  time-reversal  detection  measure  increases 
significantly. 


•  the  detection  power,  that  is  the  probability  that  the  amplitude  is  greater  than  0.91  when 
the  target  is  present,  is  over  99%  for  values  of  L  as  low  as  L  =  50  corresponding  to  a 
strong  scatterer  occupying  5%  of  the  medium; 

•  the  detection  power  decreases  for  smaller  targets  to  88%  for  L  =  40,  and  70%  for  L  =  20, 

422.  Detecting  the  location  of  the  scatterer.  The  potential  use  of  the  refocused  pulse 
amplitude  to  determine  the  location  of  the  buried  strong  scatterer  yc  was  analysed  in  the  case 
L  —  100.  We  simulated  a  time-reversal  expanding  window  experiment  based  on  increasing 
fo  by  increments  of  0.2  s  in  the  range  1-3  s  and  recording  the  statistics  of  the  refocused 
pulse  amplitude  for  ten  realizations  of  the  reference  medium  (i.e.  L  =  0).  A  relatively 
small  number  of  realizations  were  used  to  model  potential  limitations  associated  with  data 
collection  in  practical  applications.  The  expanding  window  experiment  was  then  repeated  for 
ten  independent  realizations  of  the  medium  with  the  buried  strong  scatterer  located  at  y  =  yc. 
Refocused  pulse  amplitudes  were  recorded  with  increasing  t^  in  the  cases  yc  =  0.5,  0.6,  0.7 
(figure  5).  In  all  three  cases,  it  is  observed  that  the  refocused  pulse  amplitude  exhibits  a 
statistically  significant  jump  above  the  distribution  of  reference  pulse  amplitudes  at  a  critical 
value  of  fo  that  increases  with  y^- 

Fonnulation  of  a  strategy  for  detecting  the  location  of  the  buried  scatterer  was  then 
conducted  via  analysis  of  the  expanding  window  data  for  each  realization  of  the  medium.  For 
each  realization  of  expanding  window  amplitude  data,  a  critical  expanding  window  time  fg 
was  determined  using  the  criterion: 

=  min{/o  €  [l,3]|A(ri)^=ioo  >  {A(ri)L=o}  +  2cr{A(fi)^=o},  €  [/o,  3]]  (20) 

where  A  is  the  refocused  pulse  amplitude  and  a  is  the  standard  deviation.  For  each  realization, 
this  criterion  determines  the  smallest  recording  time  at  which  the  refocused  pulse  amplitude 
exceeds  the  reference  amplitude  by  2a  in  further  expansion  of  the  recording  window.  Note  that 
for  larger  window  size  most  of  the  signal  would  be  recorded  and  consequently  the  amplitude 
of  the  refocused  pulse  would  be  the  same  with  or  without  the  buried  object.  Distributions  of  t^ 
with  respect  to  t^  (over  ten  realizations)  for  yc  =  0.5,  0.6,  0.7  are  shown  using  histograms  in 
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(a)  =0,5,  L=100,  N=1000  (b)  =0.6,  U100,  N=1000 


(c)y^=0.7,  L=100,  N=1000 


Figure  5.  Results  of  the  expanding  window  time-reversal  experiments  for  ten  independent 
realizations  of  the  medium  with  buried  object  of  size  Z.  =  100  at  y  —  yc  for  three  increasing 
values. 


Distribution  of  t^^*  in  Detection  (L=100,  N=1000) 


Figure  6.  Statistical  performance  of  the  detection  criterion  (20)  for  ten  independent  realizations 
of  the  medium  with  buried  object  of  size  L  =  100  at  y  =  y^  for  three  increasing  values.  The  most 
commonly  occurring  values  of  correspond  to  2y(;  with  a  systematic  bias  to  the  right. 


figure  6.  For  jc  =  0.5, 0.6  and  0.7,  the  most  commonly  occurring  values  of  corresponding 
to  Ijc,  are  observed  to  be  1.2, 1.4  and  1.6,  respectively.  These  critical  times  correspond  to  the 
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first  time-reversal  recording  window  that  contains  reflections  from  the  buried  scatterer.  Using 
our  criterion  (20),  detection  of  scatterer  location  is  most  accurate  in  the  case  ==  0*5,  where 
all  values  ai*e  within  one  window  increment  of  the  conect  value  of  1.2.  Accuracy  reduces  as 
yc  increases  with  8  of  10  cases  within  one  window  increment  of  1.4  for  =  0.6,  and  6  of  10 
cases  within  one  window  increment  of  1 .6  for  jc  =  0.7 . 

In  the  current  study,  we  considered  random  media  with  layers  of  fixed  width  1/A^.  To 
assess  the  potential  effects  of  varying  layer  widths,  we  computed  refocused  pulse  amplitudes  m 
the  case  /q  =  1.0  s,  L  =  0  for  50  realizations  of  the  medium.  These  realizations  were  obtained 
by  generating  a  set  of  random  cj  values  (as  described  earlier)  and  randomly  perturbing  the  fixed 
layer  widths  1/N.  Forperturbations  of  magnitude  0%,  5%,  10%  and  15%,  statistics  (MEAN± 
SD)  of  the  refocused  pulse  amplitude  were  0.7523  dr  0.0248,  0.7519  dr  0.0251,  0.7533  dr 
0.0260  and  0.7551  dr  0.0256,  respectively.  Statistics  of  the  (unsigned)  relative  error  between 
pulse  amplitudes  in  the  fixed  and  varying  layer  width  cases  show  a  mean  close  to  1%  with  a 
standard  deviation  of  less  than  1%.  This  is  in  agreement  with  the  theoretical  analysis  which 
shows  that  only  the  unsealed  integrated  autocorrelation  of  the  medium  E{Y]{0)r](x)}  dx 
plays  a  role  in  the  limit  6^0,  the  correlation  length  of  the  medium  being  of  order  6^  (see  [7]). 

These  results  show  promise  for  use  of  a  windowed  time-reversal  technique  in  detecting  the 
location  of  a  strong  scatterer  buried  in  a  randomly  layered  medium.  Inaccuracy  in  detection 
of  location  will  likely  minimize  as  the  number  of  layers  is  increased  to  enhance  self-averaging 
over  the  realizations  of  the  random  medium.  In  future  studies,  we  will  consider  larger  scale 
simulations  (N  >  10^)  and  note  that  our  frequency  domain  numerical  approach  for  simulating 
time  reversal  is  amenable  to  parallelization. 


5.  Conclusions 

Using  numerical  simulations  we  have  demonstrated  that  time-reversal  refocusing  for  acoustic 
pulses  shows  promise  for  detection  of  a  strong  scatterer  buried  in  a  random  medium.  This  is 
in  the  regime  where  the  correlation  length  of  the  medium  is  small  but  where  the  fluctuations 
are  large,  making  the  detection  problem  extiemely  challenging.  Our  numerical  experiments 
are  performed  in  the  acoustic  ID  case  with  a  method  that  is  well  adapted  to  more  realistic 
multi-dimensional  situations. 
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